In [1]:
# basic stuff
import pandas as pd
import numpy as np
from itertools import product
import holidays
import pickle
pd.set_option('display.max_columns', 500)
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)

# time stuff
from datetime import datetime
hour = pd.DateOffset(hours=1)

# modeL stuff
from sklearn.metrics import pairwise_distances_argmin_min, r2_score
from sklearn.model_selection import TimeSeriesSplit
import xgboost as xgb
from xgboost import plot_importance, XGBRegressor

# geo stuff
import geopandas as gpd
import json

# plotting stuff
import plotly
import plotly.graph_objs as go
from plotly.subplots import make_subplots
import plotly.io as pio
from plotly.offline import iplot
import matplotlib.pyplot as plt
In [2]:
pio.templates['rockwell'] = go.layout.Template(
    layout=go.Layout(
        font=dict(
            family='Rockwell',
            size=18,
            color='#2A5674'
        ),
        title={'yanchor': 'middle'},
    )
)
pio.templates.default = "plotly_white+rockwell"
In [3]:
pink= '#e2d9e2'
blue = '#6785be'
red = '#8e2c50'
orange = '#ba6657'
beige = '#ceac94'
purple = '#421257'

 1.2 Taxi zones

С сайта TLC скачиваем shapefile c координатами зон, на которые поделен город. Здесь этот файл переформатируем в geojson с координатами долготы/широты каждой зоны. Посмотрим на него еще раз:

In [4]:
zones_gdf = gpd.read_file('NYC_TAXI_data/other_data/taxi_zones.geojson').loc[:,['zone', 'OBJECTID', 'borough','geometry']]
zones_gdf.columns = ['zone_name', 'zone_id', 'borough','geometry']
zones_gdf['zone_id'] = zones_gdf['zone_id'].astype(str)
NY_center_lat = (40.49612+40.91553)/2
NY_center_lon = (-74.25559-73.70001)/2
In [5]:
zones_gdf.sample(5)
Out[5]:
zone_name zone_id borough geometry
165 Morningside Heights 166 Manhattan POLYGON ((-73.95708 40.81430, -73.95627 40.813...
103 Governor's Island/Ellis Island/Liberty Island 104 Manhattan POLYGON ((-74.03995 40.70089, -74.03945 40.700...
251 Whitestone 252 Queens POLYGON ((-73.82050 40.80101, -73.82040 40.800...
148 Madison 149 Brooklyn POLYGON ((-73.94406 40.61199, -73.94380 40.610...
93 Fordham South 94 Bronx POLYGON ((-73.89964 40.86221, -73.89944 40.862...
In [6]:
zones_geojson = json.loads(zones_gdf.to_json())
layout = go.Layout(mapbox_style='carto-positron',
                   mapbox_zoom=9.5,
                   mapbox_center = {'lat': NY_center_lat, 'lon': NY_center_lon},
                   hoverlabel=dict(
                       bgcolor="white",
                       font_size=12,
                       font_family='Rockwell'),
                   margin={"r":0,"t":50,"l":0,"b":0},
                   title='New York taxi zones',
                   height=600, width=800
                  )

data = go.Choroplethmapbox(geojson=zones_geojson,
                           locations=zones_gdf['zone_id'],
                           z=zones_gdf['zone_id'],
                           hovertext=zones_gdf['zone_name'],
                           hovertemplate='<b>Zone name</b>: <b>%{hovertext}</b>'+
                                           '<br><b>Zone ID </b>: %{location}'+
                                           "<extra></extra>",
                           showlegend=False,
                           autocolorscale=False,
                           colorscale='twilight',
                           showscale=False,
                           marker_opacity=0.8,
                           marker_line_width=0.1
                          )
taxi_zones = go.Figure(data=data, layout=layout)
In [7]:
taxi_zones.write_html("exp/interactive_taxi_zones.html")
In [8]:
taxi_zones.show()
# при наведении на зону можно увидеть ее название и ID

2. Временные ряды

Для того, чтобы решить поставленную задачу, сырые данные необходимо агрегировать по часам и районам. Агрегированные данные будут представляют собой почасовые временные ряды с количеством поездок из каждой зоны:

In [9]:
aggregated_data = pd.read_csv('NYC_TAXI_aggregated_data/pu_agg_data.csv',
                              index_col=['pickup_datetime'], parse_dates=['pickup_datetime'])
aggregated_data.index.name = 'pickup_datetime'
aggregated_data.columns.name = 'zone_id'
In [10]:
aggregated_data.head()
Out[10]:
zone_id 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265
pickup_datetime
2016-01-01 00:00:00 0.0 0.0 0.0 153.0 0.0 0.0 112.0 1.0 0.0 0.0 0.0 9.0 131.0 7.0 0.0 0.0 36.0 2.0 0.0 0.0 1.0 0.0 0.0 110.0 47.0 1.0 0.0 3.0 0.0 0.0 0.0 0.0 33.0 0.0 1.0 30.0 63.0 0.0 3.0 30.0 148.0 99.0 429.0 0.0 99.0 0.0 1.0 775.0 27.0 438.0 1.0 14.0 0.0 5.0 0.0 2.0 0.0 0.0 0.0 0.0 34.0 3.0 0.0 0.0 48.0 21.0 0.0 731.0 3.0 6.0 2.0 1.0 0.0 120.0 158.0 3.0 0.0 2.0 1244.0 75.0 0.0 19.0 7.0 0.0 2.0 0.0 196.0 54.0 14.0 574.0 0.0 1.0 0.0 5.0 5.0 0.0 38.0 0.0 0.0 194.0 0.0 0.0 0.0 0.0 0.0 19.0 808.0 0.0 0.0 0.0 0.0 70.0 363.0 468.0 0.0 92.0 0.0 0.0 3.0 0.0 0.0 0.0 1.0 0.0 106.0 3.0 16.0 0.0 31.0 3.0 0.0 191.0 4.0 2.0 0.0 6.0 362.0 49.0 0.0 291.0 643.0 591.0 269.0 295.0 71.0 58.0 1.0 590.0 0.0 0.0 274.0 37.0 2.0 0.0 0.0 0.0 6.0 393.0 9.0 0.0 700.0 426.0 425.0 709.0 0.0 139.0 1.0 13.0 4.0 804.0 3.0 0.0 4.0 2.0 0.0 0.0 6.0 0.0 43.0 2.0 113.0 2.0 0.0 0.0 0.0 433.0 0.0 4.0 32.0 8.0 0.0 0.0 14.0 1.0 6.0 6.0 1.0 3.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 1.0 0.0 1.0 48.0 1.0 212.0 2.0 2.0 0.0 1.0 1.0 9.0 0.0 0.0 1.0 1.0 1.0 41.0 93.0 22.0 40.0 5.0 18.0 618.0 189.0 422.0 146.0 272.0 924.0 4.0 554.0 828.0 529.0 676.0 0.0 1.0 0.0 36.0 83.0 0.0 385.0 7.0 0.0 664.0 2.0 0.0 0.0 0.0 1.0 137.0 122.0 3.0 0.0 0.0 26.0 97.0 340.0 723.0 0.0 0.0
2016-01-01 01:00:00 0.0 0.0 1.0 218.0 0.0 0.0 277.0 0.0 0.0 1.0 0.0 4.0 157.0 13.0 0.0 0.0 104.0 5.0 0.0 0.0 2.0 3.0 0.0 103.0 94.0 1.0 0.0 4.0 1.0 0.0 0.0 3.0 51.0 2.0 3.0 72.0 127.0 0.0 6.0 74.0 257.0 143.0 260.0 0.0 74.0 1.0 5.0 716.0 99.0 466.0 1.0 34.0 0.0 11.0 1.0 8.0 0.0 0.0 0.0 5.0 88.0 23.0 0.0 0.0 48.0 27.0 0.0 679.0 8.0 11.0 1.0 3.0 0.0 253.0 271.0 8.0 3.0 3.0 1133.0 208.0 0.0 37.0 20.0 0.0 5.0 0.0 262.0 112.0 21.0 396.0 1.0 7.0 1.0 3.0 30.0 0.0 81.0 0.0 0.0 180.0 0.0 1.0 0.0 0.0 0.0 42.0 1068.0 1.0 0.0 0.0 0.0 137.0 261.0 296.0 0.0 133.0 0.0 0.0 6.0 1.0 0.0 0.0 2.0 0.0 100.0 6.0 26.0 0.0 79.0 2.0 0.0 84.0 8.0 5.0 1.0 8.0 520.0 4.0 0.0 338.0 850.0 577.0 275.0 276.0 111.0 64.0 6.0 697.0 1.0 0.0 331.0 73.0 4.0 0.0 0.0 0.0 7.0 294.0 10.0 1.0 579.0 658.0 386.0 637.0 0.0 184.0 2.0 25.0 13.0 1094.0 1.0 0.0 10.0 6.0 0.0 0.0 6.0 0.0 91.0 0.0 241.0 1.0 0.0 0.0 2.0 556.0 0.0 26.0 72.0 9.0 0.0 0.0 13.0 0.0 15.0 7.0 4.0 20.0 0.0 0.0 0.0 6.0 0.0 0.0 0.0 1.0 0.0 0.0 60.0 0.0 159.0 1.0 4.0 0.0 0.0 3.0 8.0 0.0 0.0 13.0 0.0 0.0 106.0 163.0 48.0 100.0 4.0 46.0 887.0 156.0 502.0 260.0 342.0 532.0 6.0 672.0 819.0 625.0 795.0 0.0 3.0 0.0 65.0 158.0 0.0 323.0 22.0 2.0 419.0 2.0 0.0 0.0 0.0 1.0 225.0 232.0 8.0 2.0 1.0 52.0 124.0 354.0 780.0 0.0 0.0
2016-01-01 02:00:00 0.0 0.0 1.0 196.0 0.0 0.0 382.0 0.0 0.0 1.0 1.0 4.0 100.0 19.0 0.0 0.0 108.0 9.0 0.0 2.0 1.0 2.0 0.0 77.0 96.0 0.0 0.0 4.0 2.0 0.0 0.0 0.0 70.0 4.0 5.0 149.0 202.0 0.0 0.0 56.0 205.0 169.0 195.0 0.0 89.0 0.0 5.0 727.0 118.0 460.0 3.0 22.0 1.0 7.0 1.0 10.0 1.0 0.0 0.0 3.0 132.0 22.0 1.0 0.0 60.0 30.0 1.0 527.0 16.0 16.0 5.0 2.0 0.0 206.0 166.0 3.0 2.0 3.0 1025.0 230.0 0.0 59.0 26.0 0.0 8.0 0.0 253.0 88.0 24.0 341.0 0.0 5.0 1.0 1.0 32.0 0.0 101.0 0.0 0.0 139.0 0.0 0.0 0.0 0.0 0.0 48.0 929.0 0.0 0.0 0.0 1.0 211.0 189.0 238.0 0.0 153.0 0.0 0.0 10.0 0.0 1.0 0.0 2.0 0.0 102.0 3.0 32.0 0.0 92.0 5.0 0.0 83.0 12.0 4.0 0.0 4.0 487.0 4.0 0.0 239.0 628.0 520.0 230.0 295.0 108.0 61.0 4.0 704.0 4.0 0.0 225.0 77.0 1.0 0.0 0.0 0.0 9.0 187.0 14.0 4.0 610.0 545.0 386.0 503.0 2.0 127.0 9.0 24.0 7.0 1139.0 2.0 0.0 13.0 0.0 0.0 0.0 8.0 0.0 90.0 1.0 259.0 2.0 0.0 0.0 2.0 448.0 0.0 25.0 96.0 4.0 0.0 0.0 16.0 0.0 20.0 15.0 4.0 21.0 0.0 2.0 0.0 6.0 0.0 0.0 0.0 0.0 0.0 0.0 73.0 0.0 123.0 3.0 3.0 0.0 0.0 3.0 15.0 0.0 0.0 6.0 0.0 0.0 117.0 153.0 67.0 132.0 6.0 60.0 742.0 170.0 485.0 257.0 315.0 313.0 9.0 395.0 426.0 383.0 671.0 0.0 10.0 3.0 53.0 121.0 0.0 222.0 23.0 1.0 347.0 3.0 0.0 3.0 0.0 0.0 259.0 257.0 9.0 6.0 0.0 69.0 103.0 233.0 774.0 0.0 0.0
2016-01-01 03:00:00 0.0 0.0 0.0 199.0 0.0 1.0 286.0 0.0 0.0 0.0 2.0 5.0 44.0 21.0 0.0 2.0 95.0 2.0 0.0 1.0 3.0 2.0 1.0 66.0 65.0 3.0 0.0 0.0 1.0 0.0 0.0 1.0 50.0 1.0 8.0 120.0 180.0 0.0 2.0 50.0 210.0 144.0 73.0 0.0 67.0 0.0 0.0 959.0 92.0 463.0 2.0 19.0 0.0 7.0 0.0 8.0 1.0 0.0 1.0 1.0 94.0 5.0 1.0 0.0 37.0 10.0 2.0 504.0 5.0 16.0 1.0 6.0 0.0 126.0 129.0 6.0 4.0 1.0 1149.0 210.0 3.0 59.0 27.0 0.0 4.0 0.0 160.0 43.0 11.0 338.0 6.0 6.0 2.0 2.0 21.0 0.0 78.0 0.0 0.0 135.0 0.0 4.0 0.0 0.0 0.0 36.0 745.0 1.0 0.0 0.0 1.0 131.0 174.0 233.0 0.0 122.0 1.0 0.0 5.0 0.0 0.0 1.0 3.0 0.0 88.0 2.0 29.0 0.0 89.0 2.0 0.0 33.0 7.0 2.0 3.0 5.0 322.0 9.0 0.0 98.0 474.0 325.0 219.0 348.0 77.0 60.0 1.0 734.0 1.0 0.0 117.0 61.0 2.0 0.0 0.0 0.0 6.0 157.0 15.0 2.0 447.0 369.0 402.0 727.0 0.0 81.0 9.0 35.0 5.0 698.0 2.0 0.0 16.0 3.0 0.0 0.0 9.0 1.0 77.0 1.0 203.0 1.0 1.0 0.0 4.0 375.0 0.0 19.0 66.0 8.0 0.0 0.0 13.0 0.0 9.0 6.0 5.0 26.0 0.0 3.0 0.0 3.0 1.0 0.0 1.0 0.0 0.0 0.0 43.0 1.0 150.0 5.0 1.0 1.0 1.0 5.0 21.0 0.0 1.0 5.0 0.0 0.0 99.0 96.0 65.0 116.0 6.0 81.0 370.0 319.0 410.0 173.0 178.0 475.0 9.0 246.0 211.0 262.0 407.0 0.0 5.0 1.0 52.0 114.0 0.0 248.0 18.0 5.0 307.0 4.0 0.0 0.0 0.0 1.0 273.0 257.0 5.0 5.0 1.0 80.0 68.0 114.0 491.0 0.0 0.0
2016-01-01 04:00:00 0.0 0.0 0.0 135.0 0.0 0.0 259.0 0.0 1.0 0.0 0.0 1.0 23.0 13.0 1.0 0.0 59.0 5.0 0.0 1.0 2.0 2.0 0.0 56.0 28.0 0.0 0.0 2.0 3.0 0.0 1.0 1.0 17.0 0.0 3.0 108.0 123.0 0.0 0.0 14.0 129.0 99.0 58.0 0.0 53.0 0.0 2.0 785.0 52.0 301.0 0.0 7.0 1.0 2.0 1.0 5.0 0.0 3.0 0.0 4.0 63.0 7.0 0.0 0.0 37.0 3.0 0.0 365.0 3.0 9.0 3.0 4.0 1.0 131.0 98.0 6.0 2.0 3.0 1105.0 184.0 0.0 52.0 22.0 0.0 1.0 0.0 91.0 33.0 17.0 286.0 4.0 6.0 1.0 2.0 11.0 0.0 49.0 1.0 0.0 154.0 0.0 3.0 0.0 0.0 0.0 14.0 478.0 2.0 0.0 0.0 2.0 62.0 166.0 280.0 0.0 105.0 0.0 0.0 12.0 0.0 0.0 0.0 3.0 0.0 65.0 2.0 26.0 0.0 125.0 4.0 1.0 33.0 8.0 2.0 1.0 8.0 192.0 2.0 0.0 76.0 284.0 178.0 134.0 246.0 58.0 57.0 1.0 683.0 2.0 2.0 62.0 49.0 1.0 0.0 1.0 0.0 10.0 108.0 10.0 4.0 235.0 227.0 323.0 481.0 1.0 43.0 2.0 25.0 8.0 340.0 2.0 0.0 19.0 3.0 0.0 0.0 6.0 0.0 75.0 2.0 80.0 3.0 0.0 0.0 2.0 286.0 0.0 22.0 23.0 5.0 0.0 1.0 18.0 0.0 2.0 8.0 2.0 26.0 0.0 0.0 1.0 3.0 0.0 0.0 0.0 0.0 1.0 4.0 32.0 0.0 114.0 1.0 4.0 0.0 0.0 1.0 11.0 0.0 0.0 5.0 0.0 0.0 73.0 53.0 43.0 105.0 6.0 56.0 181.0 321.0 247.0 84.0 113.0 332.0 13.0 117.0 126.0 152.0 196.0 1.0 1.0 3.0 52.0 112.0 0.0 166.0 17.0 3.0 336.0 4.0 0.0 0.0 0.0 3.0 184.0 173.0 3.0 5.0 1.0 58.0 44.0 70.0 306.0 0.0 0.0

Нарисуем несколько временных рядов, чтобы посмотреть, как они вообще выглядят:

In [11]:
timerows_to_plot = ['7', '161', '261']
timerows_colors = {timerows_to_plot[0]: purple,
                   timerows_to_plot[1]: blue,
                   timerows_to_plot[2]: orange
                  }

interactive_some_timerows = go.Figure()
for zone in timerows_to_plot:
    df = aggregated_data[[zone]].reset_index()
    x = df['pickup_datetime']
    y = df[str(zone)]
    interactive_some_timerows.add_trace(go.Scatter(x=x, y=y,
                                                   mode='lines',
                                                   name=str(zone),
                                                   line=dict(width=2,
                                                             color=timerows_colors[zone])))
interactive_some_timerows = interactive_some_timerows.update_xaxes(rangeslider_visible=True)
interactive_some_timerows = interactive_some_timerows.update_layout(title='Timerows for selected taxi zones',
                                                                    autosize=True)
In [12]:
interactive_some_timerows.write_html("exp/interactive_some_timerows.html")
In [13]:
interactive_some_timerows.show()
# можно выбрать интересующий период на шкале снизу, а сверху все будет в увеличенном виде!

Сильно выражена суточная сезонность, а так же недельная и годовая. Ряды разные.

Обучение и тест
Нужно обозначить дату, разделяющую обучение и тест. Обучаться будем до 2018/06/30 23:00, а предсказывать будем на июль 2018 года.

In [14]:
split_date = datetime(2018,6,30,23) # время, разделяющее train и test

2.1 Отбор рядов и визуализация спроса на такси на карте

Не все ряды одинаково полезны, то есть не все зоны одинаково востребованы. Прогнозировать спрос имеет смысл в популярных зоных, чтобы понимать, сколько в зоне должно быть доступных машин в каждый час.

Посчитаем среднее число поездок из каждой зоны в час и раскрасим регионы в соответствии с полученными значениями:

In [15]:
mean_counts = pd.DataFrame(aggregated_data.loc[:split_date].mean(), columns=['mean_count']).reset_index()

zones_gdf_with_means = zones_gdf.merge(mean_counts, how='left', on='zone_id')#.fillna(0)
In [16]:
zones_gdf_with_means.head()
Out[16]:
zone_name zone_id borough geometry mean_count
0 Newark Airport 1 EWR POLYGON ((-74.18445 40.69500, -74.18449 40.695... 0.325567
1 Jamaica Bay 2 Queens MULTIPOLYGON (((-73.82338 40.63899, -73.82277 ... 0.008406
2 Allerton/Pelham Gardens 3 Bronx POLYGON ((-73.84793 40.87134, -73.84725 40.870... 0.046372
3 Alphabet City 4 Manhattan POLYGON ((-73.97177 40.72582, -73.97179 40.725... 36.748218
4 Arden Heights 5 Staten Island POLYGON ((-74.17422 40.56257, -74.17349 40.562... 0.001553
In [17]:
zones_geojson = json.loads(zones_gdf_with_means.to_json())
layout = go.Layout(
    mapbox_style='carto-positron',
    mapbox_zoom=9,
    mapbox_center = {'lat': NY_center_lat, 'lon': NY_center_lon},
    hoverlabel=dict(
        bgcolor="white",
        font_size=12,
        font_family='Rockwell'),
    margin={"r":0,"t":50,"l":0,"b":0},
    title='New York taxi zones colored by mean trip count',
    height=600, width=600
)
data = go.Choroplethmapbox(geojson=zones_geojson,
                           locations=zones_gdf_with_means['zone_id'],
                           z=zones_gdf_with_means['mean_count'],
                           hovertext=zones_gdf_with_means['zone_name'],
                           hovertemplate='<b>Zone name</b>: <b>%{hovertext}</b>'+
                                                    '<br><b>Zone ID </b>: %{location}'+
                                                    '<br><b>Mean trip count </b>: %{z}<br>'+
                                                    "<extra></extra>",
                           showlegend=False,
                           autocolorscale=False,
                           colorscale='BuPu',
                           marker_opacity=0.8, marker_line_width=0.1
                          )

mean_colored_zones = go.Figure(data=data, layout=layout)
In [18]:
mean_colored_zones.write_html("exp/mean_colored_zones.html")
In [19]:
mean_colored_zones.show()

Работать будем с зонами, из которых среднее число поездок в час не меньше 5, таких зон получется 83:

In [20]:
popular_zones = list(aggregated_data.loc[:, aggregated_data.mean() >= 5].columns)
print(len(popular_zones))
83
In [21]:
R = len(popular_zones)
In [22]:
agg_data_pop = aggregated_data[popular_zones].copy()
agg_data_pop.to_csv('NYC_TAXI_aggregated_data/pu_agg_data_pop.csv', index=True)
agg_data_pop.head()
Out[22]:
zone_id 4 7 12 13 24 25 33 40 41 42 43 45 48 49 50 52 65 66 68 74 75 79 80 87 88 90 97 100 107 112 113 114 116 125 129 132 137 138 140 141 142 143 144 145 146 148 151 152 158 161 162 163 164 166 170 179 181 186 209 211 223 224 226 229 230 231 232 233 234 236 237 238 239 244 246 249 255 256 260 261 262 263 264
pickup_datetime
2016-01-01 00:00:00 153.0 112.0 9.0 131.0 110.0 47.0 33.0 30.0 148.0 99.0 429.0 99.0 775.0 27.0 438.0 14.0 48.0 21.0 731.0 120.0 158.0 1244.0 75.0 196.0 54.0 574.0 38.0 194.0 808.0 70.0 363.0 468.0 92.0 106.0 31.0 191.0 362.0 49.0 291.0 643.0 591.0 269.0 295.0 71.0 58.0 590.0 274.0 37.0 393.0 700.0 426.0 425.0 709.0 139.0 804.0 43.0 113.0 433.0 48.0 212.0 41.0 93.0 40.0 618.0 189.0 422.0 146.0 272.0 924.0 554.0 828.0 529.0 676.0 83.0 385.0 664.0 137.0 122.0 26.0 97.0 340.0 723.0 0.0
2016-01-01 01:00:00 218.0 277.0 4.0 157.0 103.0 94.0 51.0 74.0 257.0 143.0 260.0 74.0 716.0 99.0 466.0 34.0 48.0 27.0 679.0 253.0 271.0 1133.0 208.0 262.0 112.0 396.0 81.0 180.0 1068.0 137.0 261.0 296.0 133.0 100.0 79.0 84.0 520.0 4.0 338.0 850.0 577.0 275.0 276.0 111.0 64.0 697.0 331.0 73.0 294.0 579.0 658.0 386.0 637.0 184.0 1094.0 91.0 241.0 556.0 60.0 159.0 106.0 163.0 100.0 887.0 156.0 502.0 260.0 342.0 532.0 672.0 819.0 625.0 795.0 158.0 323.0 419.0 225.0 232.0 52.0 124.0 354.0 780.0 0.0
2016-01-01 02:00:00 196.0 382.0 4.0 100.0 77.0 96.0 70.0 56.0 205.0 169.0 195.0 89.0 727.0 118.0 460.0 22.0 60.0 30.0 527.0 206.0 166.0 1025.0 230.0 253.0 88.0 341.0 101.0 139.0 929.0 211.0 189.0 238.0 153.0 102.0 92.0 83.0 487.0 4.0 239.0 628.0 520.0 230.0 295.0 108.0 61.0 704.0 225.0 77.0 187.0 610.0 545.0 386.0 503.0 127.0 1139.0 90.0 259.0 448.0 73.0 123.0 117.0 153.0 132.0 742.0 170.0 485.0 257.0 315.0 313.0 395.0 426.0 383.0 671.0 121.0 222.0 347.0 259.0 257.0 69.0 103.0 233.0 774.0 0.0
2016-01-01 03:00:00 199.0 286.0 5.0 44.0 66.0 65.0 50.0 50.0 210.0 144.0 73.0 67.0 959.0 92.0 463.0 19.0 37.0 10.0 504.0 126.0 129.0 1149.0 210.0 160.0 43.0 338.0 78.0 135.0 745.0 131.0 174.0 233.0 122.0 88.0 89.0 33.0 322.0 9.0 98.0 474.0 325.0 219.0 348.0 77.0 60.0 734.0 117.0 61.0 157.0 447.0 369.0 402.0 727.0 81.0 698.0 77.0 203.0 375.0 43.0 150.0 99.0 96.0 116.0 370.0 319.0 410.0 173.0 178.0 475.0 246.0 211.0 262.0 407.0 114.0 248.0 307.0 273.0 257.0 80.0 68.0 114.0 491.0 0.0
2016-01-01 04:00:00 135.0 259.0 1.0 23.0 56.0 28.0 17.0 14.0 129.0 99.0 58.0 53.0 785.0 52.0 301.0 7.0 37.0 3.0 365.0 131.0 98.0 1105.0 184.0 91.0 33.0 286.0 49.0 154.0 478.0 62.0 166.0 280.0 105.0 65.0 125.0 33.0 192.0 2.0 76.0 284.0 178.0 134.0 246.0 58.0 57.0 683.0 62.0 49.0 108.0 235.0 227.0 323.0 481.0 43.0 340.0 75.0 80.0 286.0 32.0 114.0 73.0 53.0 105.0 181.0 321.0 247.0 84.0 113.0 332.0 117.0 126.0 152.0 196.0 112.0 166.0 336.0 184.0 173.0 58.0 44.0 70.0 306.0 0.0
In [23]:
del aggregated_data
In [24]:
popmean_counts = pd.DataFrame(agg_data_pop.mean(), columns=['mean_count']).reset_index()

popzones_with_means = zones_gdf.merge(popmean_counts, how='left', on='zone_id')
popzones_with_means.head()
Out[24]:
zone_name zone_id borough geometry mean_count
0 Newark Airport 1 EWR POLYGON ((-74.18445 40.69500, -74.18449 40.695... NaN
1 Jamaica Bay 2 Queens MULTIPOLYGON (((-73.82338 40.63899, -73.82277 ... NaN
2 Allerton/Pelham Gardens 3 Bronx POLYGON ((-73.84793 40.87134, -73.84725 40.870... NaN
3 Alphabet City 4 Manhattan POLYGON ((-73.97177 40.72582, -73.97179 40.725... 34.129486
4 Arden Heights 5 Staten Island POLYGON ((-74.17422 40.56257, -74.17349 40.562... NaN
In [25]:
popzones_geojson = json.loads(popzones_with_means.to_json())

layout = go.Layout(mapbox_style='carto-positron',
                                          mapbox_zoom=9,
                                          mapbox_center = {'lat': NY_center_lat, 'lon': NY_center_lon},
                                          hoverlabel=dict(
                                              bgcolor="white",
                                              font_size=12,
                                              font_family='Rockwell'),
                                          margin={"r":0,"t":50,"l":0,"b":0},
                                          title='Most popular New York taxi zones<br>colored by mean trip count',
                                                             height=600, width=600)
data =  go.Choroplethmapbox(geojson=popzones_geojson,
                                    locations=popzones_with_means['zone_id'],
                                    z=popzones_with_means['mean_count'],
                                    hovertext=popzones_with_means['zone_name'],
                                    hovertemplate='<b>Zone name</b>: <b>%{hovertext}</b>'+
                                                    '<br><b>Zone ID </b>: %{location}'+
                                                    '<br><b>Mean trip count </b>: %{z}<br>'+
                                                    "<extra></extra>",
                                    showlegend=False,
                                    autocolorscale=False,
                                    colorscale='BuPu',
                                    marker_opacity=0.8, marker_line_width=0.1
                                   )
mean_colored_pop_zones = go.Figure(data=data, layout=layout)
In [26]:
mean_colored_pop_zones.show()

2.2. Кластеризация рядов

Так как подбирать параметры для отдельных моделей каждой географической зоны относительно долго, то кластеризуем ряды и подбирать параметры модели будем уже для каждого кластера отдельно.
Кластеризацию делаю методом KMeans.

In [27]:
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn.decomposition import PCA
In [28]:
# стандартизуем временные ряды
scaler = StandardScaler()
data_for_clusterization = pd.DataFrame(scaler.fit_transform(agg_data_pop.loc[:split_date]),
                                   columns=agg_data_pop.loc[:split_date].columns,
                                   index=agg_data_pop.loc[:split_date].index).T
In [29]:
clusterization = make_subplots(rows=2, cols=3, subplot_titles=('n_clusters = 2',
                                                    'n_clusters = 3',
                                                    'n_clusters = 4',
                                                    'n_clusters = 5',
                                                    'n_clusters = 6',
                                                    'n_clusters = 7',
                                                   ))
for n_clusters, (row, col) in zip(range(2,8), product([1,2],[1,2,3])):
    # кластеризуем
    clusterer = KMeans(n_clusters=n_clusters, random_state=21)
    labels = clusterer.fit_predict(data_for_clusterization)
    centers = clusterer.cluster_centers_
    # уменьшим размерность
    reductor = PCA(2, random_state=21)
    dots = reductor.fit_transform(data_for_clusterization)
    dots = pd.DataFrame(dots, columns=['x','y'])
    dots['labels'] = labels.astype(float)
    # добавим 2D-график
    clusterization.add_trace(go.Scatter(x=dots['x'], y=dots['y'], mode='markers',
                                        marker=dict(
                                            size=10,
                                            line_width=0.5,
                                            color=dots['labels'],
                                            colorscale='Viridis',
                                            showscale=False,
                                        ),
                                        text=dots['labels'],
                                        hovertemplate='<b>Cluster #%{text}</b>'+'<extra></extra>'
                                       ),
                             row=row, col=col)
    clusterization.update_xaxes(title_text='silhouette_score {:.2f}'.format(silhouette_score(data_for_clusterization, labels)),
                     title_font=dict(size=14, family='Rockwell'),
                     row=row, col=col)

clusterization = clusterization.update_layout(height=600, width=800, title_text='Clusterization of popular zones', showlegend=False)
clusterization = clusterization.update_xaxes(showticklabels=False)
clusterization = clusterization.update_yaxes(showticklabels=False)
In [30]:
clusterization.write_html("exp/clusterization.html")
In [31]:
clusterization.show()
# кластеризация стандартизованных временных рядов методом KMeans, подсчет коэффициента силуэта

Делить будем на 3 кластера.

In [32]:
n_clusters = 3
clusterer = KMeans(n_clusters=n_clusters, random_state=21)
labels = clusterer.fit_predict(data_for_clusterization)
centers = clusterer.cluster_centers_
In [33]:
zone_cluster_labels = {zone:label for zone, label in zip(popular_zones, labels)}

Найдем ряды, ближайшие к центрам кластеров:

In [34]:
# центры кластеров
closest, _ = pairwise_distances_argmin_min(centers, data_for_clusterization)
centroid_timerows = agg_data_pop.iloc[:,closest].copy()
centroid_names = centroid_timerows.columns.values # типичные ряды для каждого кластера
with open('XGBoost/centroid_names.pickle','wb') as file:
    pickle.dump(centroid_names, file)
In [35]:
del data_for_clusterization
In [36]:
print('Ближайшие к центрам кластеров ряды: {}'.format(centroid_names))
Ближайшие к центрам кластеров ряды: ['170' '79' '42']
In [37]:
timerows_colors = {centroid_names[0]: purple,
                   centroid_names[1]: blue,
                   centroid_names[2]: orange,
                }

centroid_timerows = go.Figure()
for zone in centroid_names:
    x = agg_data_pop[[zone]].index
    y = agg_data_pop[[zone]][str(zone)]
    centroid_timerows.add_trace(go.Scatter(x=x, y=y, mode='lines',
                                           name=str(zone),
                                           line=dict(width=2, color=timerows_colors[zone])))
centroid_timerows = centroid_timerows.update_xaxes(rangeslider_visible=True)
centroid_timerows = centroid_timerows.update_layout(title='Timerows for "centroid" taxi zones',
                                                    autosize=True)
In [38]:
centroid_timerows.write_html("exp/centroid_timerows.html")
In [39]:
centroid_timerows.show()

3. Подготовка выборок для разных задач прогнозирования

Для каждой из шести задач прогнозирования $\hat{y}_{T+i|T}, i=1,\dots,6$ будем формировать выборки, в которых отклик = $y_{T+i}$ при всевозможных значениях $T$ (то есть двигаем значения ряда "назад").

In [40]:
shifted_data = {}
for shift in range(1,7):
    shifted_data[shift] = agg_data_pop.shift(-shift).dropna()

4. Построение признаковых описаний рядов

Теперь к временному ряду каждого региона добавим признаки. Признаки берутся из:

  1. Значений времени ряда:
    • год, месяц, день месяца, день недели, час
    • праздник/обычный день
    • сезонные гармоники
  2. Признаки, зависящие от значений ряда (авторегрессионные):

    • cреднее количество поездок за предыдущие 12/24/48/168 часов
    • суммарное количество поездок за предыдущие 12/24/48/168 часов
    • среднее число поездок в такой час/день/день недели/выходной
    • количество поездок из рассматриваемого района в моменты времени $y_{T-1}, ... , y_{T-K}$
    • количество поездок из рассматриваемого района в моменты времени $y_{T-24}, y_{T-48}, ..., y_{T-24*K_d}$
  3. Статистической информации о регионе из сырых данных (эти признаки напрямую не зависят от значений целевой переменной):

    • средняя длительность поездок
    • среднее количество пассажиров
    • среднее расстояние по счётчику
    • доли географических зон, в которые совершаются поездки
    • доли поездок, совершаемых по тарифам каждого из типов
    • доли способов оплаты поездок
    • средняя стоимость поездок
    • доли провайдеров данных

4.1. Метаданные для создания признаков

Загружаем размеченные по регионам записи о поездках - marked_data за каждый месяц. Я не могу хранить такие данные в памяти, поэтому надо подумать здесь о БД!

marked_data = pd.read_csv('NYC_TAXI_marked_data/marked_data_2016_2019.csv', # вот это будет очень большой файл index_col=['pickup_datetime'], parse_dates=['pickup_datetime']) marked_data.index.name = 'pickup_datetime' marked_data['PU_zone_id'] = marked_data['PU_zone_id'].astype(str) marked_data.head()

4.2. Функция, для добавления признаков в датафрейм с временным рядом

Признаков делаем много, потом проведем отбор.

In [41]:
def add_regressors(timeseries, K_s, K_h, K_d, meta_data=None):
    timerow = timeseries.copy() # сюда будем записывать регрессоры
    timerow['datetime'] = timerow.index # вспомогательный столбец
    
    ### 1.1 Признаки, зависящие от datetime
    # год, месяц, день месяца, день недели, час
    timerow['year'] = timerow['datetime'].dt.year
    timerow['month'] = timerow['datetime'].dt.month
    timerow['quarter'] = timerow['datetime'].dt.quarter
    timerow['day'] = timerow['datetime'].dt.day
    timerow['hour'] = timerow['datetime'].dt.hour
    timerow['weekday'] = timerow['datetime'].dt.dayofweek
    timerow['dayofyear'] = timerow['datetime'].dt.dayofyear
    timerow['weekofyear'] = timerow['datetime'].dt.week
    timerow['weekend'] = (timerow['datetime'].dt.weekday>=5).astype(int)

    # праздник/обычный день
    us_holidays = holidays.UnitedStates() # данные об американских праздниках
    timerow['holiday'] = np.array([date in us_holidays for date in timerow['datetime']]).astype(int)
    
    # сезонные гармоники
    T=len(timerow.index)
    for i in range(1,K_s+1):
        timerow['dsin_'+str(i)] = np.sin(np.linspace(1,T,T)*((2*np.pi*i)/24.))
        timerow['dcos_'+str(i)] = np.cos(np.linspace(1,T,T)*((2*np.pi*i)/24.))
        timerow['wsin_'+str(i)] = np.sin(np.linspace(1,T,T)*((2*np.pi*i)/168.))
        timerow['wcos_'+str(i)] = np.cos(np.linspace(1,T,T)*((2*np.pi*i)/168.))
        timerow['ysin_'+str(i)] = np.sin(np.linspace(1,T,T)*((2*np.pi*i)/8766.))
        timerow['ycos_'+str(i)] = np.cos(np.linspace(1,T,T)*((2*np.pi*i)/8766.))
        
    ### 1.2 Признаки, зависящие от предыдущих значений ряда
    #cреднее количество поездок за предыдущие 12/24/48/168 часов
    timerow['prev_mean_12'] = timerow['real_counts'].rolling(window=12).mean()
    timerow['prev_mean_24'] = timerow['real_counts'].rolling(window=24).mean()
    timerow['prev_mean_48'] = timerow['real_counts'].rolling(window=48).mean()
    timerow['prev_mean_168'] = timerow['real_counts'].rolling(window=168).mean()
    # суммарное количество поездок за предыдущие 12/24/48/168 часов
    timerow['prev_sum_12'] = timerow['real_counts'].rolling(window=12).sum()
    timerow['prev_sum_24'] = timerow['real_counts'].rolling(window=24).sum()
    timerow['prev_sum_48'] = timerow['real_counts'].rolling(window=48).sum()
    timerow['prev_sum_168'] = timerow['real_counts'].rolling(window=168).sum()
    # среднее число поездок в такой час/день/день недели/выходной
    timerow['hour_mean'] = timerow.groupby(by=['hour'])['real_counts'].transform('mean')
    timerow['day_mean'] = timerow.groupby(by=['day'])['real_counts'].transform('mean')
    timerow['dayofweek_mean'] = timerow.groupby(by=['weekday'])['real_counts'].transform('mean')
    timerow['weekend_mean'] = timerow.groupby(by=['weekend'])['real_counts'].transform('mean')
    
    # количество поездок из рассматриваемого района в моменты времени $y_{T-1}, ... , y_{T-K}$
    for i in range(1,K_h+1): 
        timerow['prev_count_'+str(i)] = timerow['real_counts'].shift(i)
    
    # количество поездок из рассматриваемого района в моменты времени $y_{T-24}, y_{T-48}, ..., y_{T-24*K_d}$
    for i in range(1,K_d+1):
        timerow['prev_count_'+str(i*24)] = timerow['real_counts'].shift(i*24)
        
    ### 1.3 Признаки, зависящие от метаданных географической зоны
    # идентификатор зоны
    zone = timeseries.columns.name
    timerow['zone'] = int(zone)
    
    #a = meta_data.groupby(by=['pickup_datetime']).agg(np.mean)
    #timerow['mean_passenger_count'] = a['passenger_count']# среднее количество пассажиров
    #timerow['mean_trip_distance'] = a['trip_distance']# среднее расстояние по счётчику
    #timerow['mean_trip_duration'] = a['trip_duration']# средняя длительность поездок
    #timerow['mean_total_cost'] = a['total_cost']# средняя стоимость поездок
    
    # доли способов оплаты поездок
    #for x in np.unique(meta_data['payment_type']):
    #    timerow['payment_type_{}_ratio'.format(int(x))] = meta_data[meta_data['payment_type'] == x]['payment_type'].groupby(by=['pickup_datetime']).count()\
    #                                        / timerow['real_counts']
 
    # доли провайдеров данных
    #for x in np.unique(meta_data['vendorID']):
    #    timerow['vendor_{}_ratio'.format(int(x))] = meta_data[meta_data['vendorID'] == x]['vendorID'].groupby(by=['pickup_datetime']).count()\
    #                                        / timerow['real_counts']
    
    timerow.drop('datetime', axis=1, inplace=True)
    
    timerow.fillna(0, inplace=True)
    
    return timerow
In [42]:
#по умолчанию
K_s, K_h, K_d = (30,14,7) #сделаем много авторегрессионных признаков, а потом устроим отбор
In [43]:
zone = centroid_names[0]
timeseries = agg_data_pop.loc[:, [zone]]
timeseries.columns = ['real_counts']
timeseries.columns.name = zone
#zone_data = marked_data[marked_data['PU_zone_id'].astype(str) == zone]

timeseries_with_features = add_regressors(timeseries, K_s, K_h, K_d)
timeseries_with_features.columns.name = zone

После применения функции добавления признаков к ряду получается такой датафрейм:

In [44]:
timeseries_with_features.sample(5)
Out[44]:
170 real_counts year month quarter day hour weekday dayofyear weekofyear weekend holiday dsin_1 dcos_1 wsin_1 wcos_1 ysin_1 ycos_1 dsin_2 dcos_2 wsin_2 wcos_2 ysin_2 ycos_2 dsin_3 dcos_3 wsin_3 wcos_3 ysin_3 ycos_3 dsin_4 dcos_4 wsin_4 wcos_4 ysin_4 ycos_4 dsin_5 dcos_5 wsin_5 wcos_5 ysin_5 ycos_5 dsin_6 dcos_6 wsin_6 wcos_6 ysin_6 ycos_6 dsin_7 dcos_7 wsin_7 wcos_7 ysin_7 ycos_7 dsin_8 dcos_8 wsin_8 wcos_8 ysin_8 ycos_8 dsin_9 dcos_9 wsin_9 wcos_9 ysin_9 ycos_9 dsin_10 dcos_10 wsin_10 wcos_10 ysin_10 ycos_10 dsin_11 dcos_11 wsin_11 wcos_11 ysin_11 ycos_11 dsin_12 dcos_12 wsin_12 wcos_12 ysin_12 ycos_12 dsin_13 dcos_13 wsin_13 wcos_13 ysin_13 ycos_13 dsin_14 dcos_14 wsin_14 wcos_14 ysin_14 ycos_14 dsin_15 dcos_15 wsin_15 wcos_15 ysin_15 ycos_15 dsin_16 dcos_16 wsin_16 wcos_16 ysin_16 ycos_16 dsin_17 dcos_17 wsin_17 wcos_17 ysin_17 ycos_17 dsin_18 dcos_18 wsin_18 wcos_18 ysin_18 ycos_18 dsin_19 dcos_19 wsin_19 wcos_19 ysin_19 ycos_19 dsin_20 dcos_20 wsin_20 wcos_20 ysin_20 ycos_20 dsin_21 dcos_21 wsin_21 wcos_21 ysin_21 ycos_21 dsin_22 dcos_22 wsin_22 wcos_22 ysin_22 ycos_22 dsin_23 dcos_23 wsin_23 wcos_23 ysin_23 ycos_23 dsin_24 dcos_24 wsin_24 wcos_24 ysin_24 ycos_24 dsin_25 dcos_25 wsin_25 wcos_25 ysin_25 ycos_25 dsin_26 dcos_26 wsin_26 wcos_26 ysin_26 ycos_26 dsin_27 dcos_27 wsin_27 wcos_27 ysin_27 ycos_27 dsin_28 dcos_28 wsin_28 wcos_28 ysin_28 ycos_28 dsin_29 dcos_29 wsin_29 wcos_29 ysin_29 ycos_29 dsin_30 dcos_30 wsin_30 wcos_30 ysin_30 ycos_30 prev_mean_12 prev_mean_24 prev_mean_48 prev_mean_168 prev_sum_12 prev_sum_24 prev_sum_48 prev_sum_168 hour_mean day_mean dayofweek_mean weekend_mean prev_count_1 prev_count_2 prev_count_3 prev_count_4 prev_count_5 prev_count_6 prev_count_7 prev_count_8 prev_count_9 prev_count_10 prev_count_11 prev_count_12 prev_count_13 prev_count_14 prev_count_24 prev_count_48 prev_count_72 prev_count_96 prev_count_120 prev_count_144 prev_count_168 zone
pickup_datetime
2016-05-30 06:00:00 88.0 2016 5 2 30 6 0 151 22 0 1 0.965926 -0.258819 0.185912 -0.982566 0.527973 -0.849261 -0.500000 -8.660254e-01 -0.365341 0.930874 -0.896774 0.442490 -0.707107 7.071068e-01 0.532032 -0.846724 0.995218 0.097683 8.660254e-01 0.5 -0.680173 0.733052 -0.793626 -0.608406 0.258819 -0.965926 0.804598 -0.593820 0.352774 0.935709 -1.000000e+00 -5.921658e-14 -0.900969 4.338837e-01 0.194431 -0.980916 0.258819 0.965926 0.965926 -0.258819 -0.683020 0.730400 8.660254e-01 -0.5 -0.997204 0.074730 0.965694 -0.259684 -0.707107 -7.071068e-01 0.993712 0.111964 -0.957232 -0.289320 -0.500000 8.660254e-01 -0.955573 -0.294755 0.660187 0.751101 0.965926 0.258819 0.884115 0.467269 -0.164111 -0.986442 1.184332e-13 -1.0 -7.818315e-01 -0.623490 -0.381441 0.924393 -0.965926 0.258819 0.652287 0.757972 0.811998 -0.583661 0.500000 8.660254e-01 -0.500000 -8.660254e-01 -0.997755 0.066968 0.707107 -7.071068e-01 0.330279 0.943883 0.882712 0.469915 -8.660254e-01 -0.5 -0.149042 -0.988831 -0.501551 -0.865128 -0.258819 0.965926 -0.037391 0.999301 -0.030816 0.999525 1.000000e+00 1.996639e-12 0.222521 -9.749279e-01 0.553893 -0.832588 -0.258819 -0.965926 -0.399892 0.916562 -0.909984 0.414644 -8.660254e-01 0.5 0.563320 -0.826239 0.991735 0.128305 0.707107 7.071068e-01 -0.707107 7.071068e-01 -0.774500 -0.632573 0.500000 -8.660254e-01 0.826239 -0.563320 0.323772 0.946135 -0.965926 -0.258819 -0.916562 0.399892 0.224567 -0.974459 -2.368663e-13 1.0 9.749279e-01 -0.222521 -0.705204 0.709005 0.965926 -0.258819 -0.999301 0.037391 0.973237 -0.229802 -0.500000 -8.660254e-01 0.988831 0.149042 -0.947862 -0.318681 -0.707107 7.071068e-01 -0.943883 -0.330279 0.636728 0.771089 8.660254e-01 0.5 8.660254e-01 0.5 -0.133635 -0.991031 0.258819 -0.965926 -0.757972 -0.652287 -0.409747 0.912199 -1.000000e+00 -5.753051e-12 0.623490 7.818315e-01 0.829598 -0.558361 236.333333 313.625000 340.041667 472.083333 2836.0 7527.0 16322.0 79310.0 225.134124 405.210859 392.693471 448.311754 56.0 51.0 67.0 142.0 204.0 255.0 337.0 405.0 441.0 349.0 441.0 478.0 552.0 392.0 89.0 129.0 334.0 383.0 358.0 432.0 367.0 170
2017-01-09 18:00:00 674.0 2017 1 1 9 18 0 9 2 0 0 -0.965926 0.258819 -0.258819 -0.965926 0.163404 0.986559 -0.500000 -8.660254e-01 0.500000 0.866025 0.322415 0.946598 0.707107 -7.071068e-01 -0.707107 -0.707107 0.472759 0.881192 8.660254e-01 0.5 0.866025 0.500000 0.610395 0.792097 -0.258819 0.965926 -0.965926 -0.258819 0.731623 0.681710 -1.000000e+00 2.290358e-13 1.000000 1.621723e-13 0.833183 0.552997 -0.258819 -0.965926 -0.965926 0.258819 0.912346 0.409420 8.660254e-01 -0.5 0.866025 -0.500000 0.966984 0.254836 0.707107 7.071068e-01 -0.707107 0.707107 0.995628 0.093402 -0.500000 8.660254e-01 0.500000 -0.866025 0.997509 -0.070543 -0.965926 -0.258819 -0.258819 0.965926 0.972575 -0.232591 -4.580715e-13 -1.0 3.243447e-13 -1.000000 0.921496 -0.388388 0.965926 -0.258819 0.258819 0.965926 0.845647 -0.533743 0.500000 8.660254e-01 -0.500000 -8.660254e-01 0.747065 -0.664751 -0.707107 7.071068e-01 0.707107 0.707107 0.628401 -0.777890 -8.660254e-01 -0.5 -0.866025 -0.500000 0.492845 -0.870117 0.258819 -0.965926 0.965926 0.258819 0.344040 -0.938955 1.000000e+00 -6.871073e-13 -1.000000 -3.176963e-14 0.185987 -0.982552 0.258819 0.965926 0.965926 -0.258819 0.022935 -0.999737 -8.660254e-01 0.5 -0.866025 0.500000 -0.140734 -0.990047 -0.707107 -7.071068e-01 0.707107 -7.071068e-01 -0.300620 -0.953744 0.500000 -8.660254e-01 -0.500000 0.866025 -0.452425 -0.891802 0.965926 0.258819 0.258819 -0.965926 -0.592068 -0.805888 9.161430e-13 1.0 -6.486893e-13 1.000000 -0.715795 -0.698310 -0.965926 0.258819 -0.258819 -0.965926 -0.820281 -0.571961 -0.500000 -8.660254e-01 0.500000 0.866025 -0.902716 -0.430236 0.707107 -7.071068e-01 -0.707107 -0.707107 -0.960885 -0.276946 8.660254e-01 0.5 8.660254e-01 0.5 -0.993224 -0.116212 -0.258819 0.965926 -0.965926 -0.258819 -0.998864 0.047647 -1.000000e+00 -6.130779e-12 1.000000 3.561143e-13 -0.977653 0.210225 580.166667 379.083333 348.666667 382.303571 6962.0 9098.0 16736.0 64227.0 636.177920 430.678241 392.693471 448.311754 599.0 466.0 501.0 536.0 570.0 583.0 560.0 508.0 701.0 697.0 567.0 304.0 117.0 53.0 352.0 393.0 743.0 675.0 650.0 600.0 442.0 170
2016-12-23 12:00:00 605.0 2016 12 4 23 12 4 358 51 0 0 -0.258819 -0.965926 0.467269 0.884115 -0.132214 0.991221 0.500000 8.660254e-01 0.826239 0.563320 -0.262106 0.965039 -0.707107 -7.071068e-01 0.993712 0.111964 -0.387397 0.921913 8.660254e-01 0.5 0.930874 -0.365341 -0.505885 0.862601 -0.965926 -0.258819 0.652287 -0.757972 -0.615492 0.788143 1.000000e+00 8.362972e-13 0.222521 -9.749279e-01 -0.714292 0.699848 -0.965926 0.258819 -0.258819 -0.965926 -0.800551 0.599265 8.660254e-01 -0.5 -0.680173 -0.733052 -0.872754 0.488160 -0.707107 7.071068e-01 -0.943883 -0.330279 -0.929634 0.368484 0.500000 -8.660254e-01 -0.988831 0.149042 -0.970192 0.242339 -0.258819 0.965926 -0.804598 0.593820 -0.993715 0.111939 1.672594e-12 -1.0 -4.338837e-01 0.900969 -0.999791 -0.020426 0.258819 0.965926 0.037391 0.999301 -0.988314 -0.152433 -0.500000 -8.660254e-01 0.500000 8.660254e-01 -0.959484 -0.281764 0.707107 7.071068e-01 0.846724 0.532032 -0.913808 -0.406147 -8.660254e-01 -0.5 0.997204 0.074730 -0.852087 -0.523400 0.965926 0.258819 0.916562 -0.399892 -0.775406 -0.631463 -1.000000e+00 -2.508891e-12 0.623490 -7.818315e-01 -0.685111 -0.728438 0.965926 -0.258819 0.185912 -0.982566 -0.582787 -0.812625 -8.660254e-01 0.5 -0.294755 -0.955573 -0.470231 -0.882543 0.707107 -7.071068e-01 -0.707107 -7.071068e-01 -0.349418 -0.936967 -0.500000 8.660254e-01 -0.955573 -0.294755 -0.222471 -0.974939 0.258819 -0.965926 -0.982566 0.185912 -0.091618 -0.995794 -3.345189e-12 1.0 -7.818315e-01 0.623490 0.040844 -0.999166 -0.258819 -0.965926 -0.399892 0.916562 0.172589 -0.984994 0.500000 8.660254e-01 0.074730 0.997204 0.301304 -0.953528 -0.707107 -7.071068e-01 0.532032 0.846724 0.424728 -0.905321 8.660254e-01 0.5 8.660254e-01 0.5 0.540696 -0.841218 -0.965926 -0.258819 0.999301 0.037391 0.647170 -0.762346 1.000000e+00 1.145744e-11 0.900969 -4.338837e-01 0.742281 -0.670089 317.916667 459.875000 472.270833 463.446429 3815.0 11037.0 22669.0 77859.0 558.466241 396.304398 463.848992 448.311754 540.0 475.0 455.0 441.0 310.0 217.0 117.0 94.0 96.0 192.0 273.0 403.0 502.0 563.0 683.0 688.0 688.0 636.0 598.0 669.0 665.0 170
2016-06-04 20:00:00 691.0 2016 6 2 4 20 5 156 22 1 0 -0.707107 0.707107 0.993712 -0.111964 0.444096 -0.895979 -1.000000 -2.984013e-13 -0.222521 -0.974928 -0.795801 0.605558 -0.707107 -7.071068e-01 -0.943883 0.330279 0.981947 -0.189155 5.968027e-13 -1.0 0.433884 0.900969 -0.963807 -0.266599 0.707107 -0.707107 0.846724 -0.532032 0.745156 0.666890 1.000000e+00 6.678303e-13 -0.623490 -7.818315e-01 -0.371481 -0.928441 0.707107 0.707107 -0.707107 0.707107 -0.079477 0.996837 -1.193605e-12 1.0 0.781831 0.623490 0.513901 -0.857849 -0.707107 7.071068e-01 0.532032 -0.846724 -0.841412 0.540394 -1.000000 1.009104e-12 -0.900969 -0.433884 0.993875 -0.110514 -0.707107 -0.707107 -0.330279 0.943883 -0.939570 -0.342357 1.335661e-12 -1.0 9.749279e-01 0.222521 0.689796 0.724004 0.707107 -0.707107 0.111964 -0.993712 -0.296516 -0.955028 1.000000 1.861436e-12 -1.000000 -2.984013e-13 -0.158452 0.987367 0.707107 7.071068e-01 0.111964 0.993712 0.580455 -0.814292 -2.387211e-12 1.0 0.974928 -0.222521 -0.881699 0.471812 -0.707107 0.707107 -0.330279 -0.943883 0.999514 -0.031174 -1.000000e+00 7.249931e-13 -0.900969 4.338837e-01 -0.909388 -0.415948 -0.707107 -0.707107 0.532032 0.846724 0.630072 0.776537 -2.018208e-12 -1.0 0.781831 -0.623490 -0.219675 -0.975573 0.707107 -7.071068e-01 -0.707107 -7.071068e-01 -0.236424 0.971650 1.000000 3.265569e-13 -0.623490 0.781831 0.643337 -0.765584 0.707107 0.707107 0.846724 0.532032 -0.916409 0.400244 -2.671321e-12 1.0 4.338837e-01 -0.900969 0.998830 0.048363 -0.707107 0.707107 -0.943883 -0.330279 -0.873453 -0.486908 -1.000000 -1.378107e-12 -0.222521 0.974928 0.566362 0.824157 -0.707107 -7.071068e-01 0.993712 0.111964 -0.141444 -0.989946 3.722871e-12 -1.0 5.968027e-13 -1.0 -0.312900 0.949786 0.707107 -0.707107 -0.993712 0.111964 0.702148 -0.712031 1.000000e+00 6.067636e-12 0.222521 9.749279e-01 -0.945320 0.326144 548.583333 491.416667 533.562500 458.458333 6583.0 11794.0 25611.0 77021.0 607.114964 403.539352 382.394374 351.717622 650.0 590.0 531.0 448.0 458.0 530.0 554.0 634.0 573.0 480.0 444.0 319.0 236.0 126.0 838.0 938.0 789.0 706.0 409.0 349.0 402.0 170
2016-10-06 21:00:00 799.0 2016 10 4 6 21 3 280 40 0 0 -0.500000 0.866025 -0.074730 0.997204 -0.994715 0.102675 -0.866025 5.000000e-01 -0.149042 0.988831 -0.204265 -0.978916 -1.000000 -8.666352e-14 -0.222521 0.974928 0.952769 -0.303695 -8.660254e-01 -0.5 -0.294755 0.955573 0.399915 0.916552 -0.500000 -0.866025 -0.365341 0.930874 -0.870647 0.491909 1.733270e-13 -1.000000e+00 -0.433884 9.009689e-01 -0.578703 -0.815539 0.500000 -0.866025 -0.500000 0.866025 0.751810 -0.659379 8.660254e-01 -0.5 -0.563320 0.826239 0.733087 0.680135 1.000000 1.169485e-12 -0.623490 0.781831 -0.601271 0.799045 0.866025 5.000000e-01 -0.680173 0.733052 -0.856558 -0.516051 0.500000 0.866025 -0.733052 0.680173 0.425377 -0.905016 -3.466541e-13 1.0 -7.818315e-01 0.623490 0.943909 0.330207 -0.500000 0.866025 -0.826239 0.563320 -0.231546 0.972824 -0.866025 5.000000e-01 -0.866025 5.000000e-01 -0.991457 -0.130437 -1.000000 -3.161802e-12 -0.900969 0.433884 0.027950 -0.999609 -8.660254e-01 -0.5 -0.930874 0.365341 0.997196 -0.074832 -0.500000 -0.866025 -0.955573 0.294755 0.176824 0.984243 2.338971e-12 -1.000000e+00 -0.974928 2.225209e-01 -0.960885 0.276946 0.500000 -0.866025 -0.988831 0.149042 -0.374141 -0.927372 8.660254e-01 -0.5 -0.997204 0.074730 0.884056 -0.467382 1.000000 -2.121839e-12 -1.000000 -8.666352e-14 0.555682 0.831395 0.866025 5.000000e-01 -0.997204 -0.074730 -0.769946 0.638109 0.500000 0.866025 -0.988831 -0.149042 -0.713790 -0.700359 -6.933082e-13 1.0 -9.749279e-01 -0.222521 0.623370 -0.781927 -0.500000 0.866025 -0.955573 -0.294755 0.841799 0.539791 -0.866025 5.000000e-01 -0.930874 -0.365341 -0.450506 0.892773 -1.000000 -3.508456e-12 -0.900969 -0.433884 -0.934311 -0.356460 -8.660254e-01 -0.5 -8.660254e-01 -0.5 0.258646 -0.965972 -0.500000 -0.866025 -0.826239 -0.563320 0.987423 0.158098 6.323603e-12 -1.000000e+00 -0.781831 -6.234898e-01 -0.055879 0.998438 658.583333 515.833333 498.979167 441.309524 7903.0 12380.0 23951.0 74140.0 557.433394 434.813657 479.684028 448.311754 855.0 754.0 722.0 503.0 480.0 534.0 620.0 609.0 746.0 596.0 685.0 794.0 931.0 660.0 682.0 600.0 506.0 386.0 487.0 726.0 794.0 170

Данные подготовили, дальше - модель. Я выбрала XGBoost (быстро, удобно, качественно).

5. Отбор признаков для XGBoost

Некоторые признаки могут быть несущественными или незначимыми для целевой переменной, и их включение в модель приводит к увеличению сложности модели, которую сложнее интерпретировать, а так же больше количество признаков увеличивает время обучения модели и уменьшает надежность предсказаний.

In [45]:
# define custom class to fix bug in xgboost 1.1.0
class myXGBR(XGBRegressor):
    @property
    def coef_(self):
        return None

Отбор признаков проводила следующим образом. Для центроида каждого кластера по кросс-валидации для временных рядов [1] оценивала качество модели, обученной на всех созданных признаках, по метрике r2_score, и фильтровала признаки по значению features_importances (XGBR.featureimportances) с порогом 0.01. Затем обучала модель на отобранных признаках и снова считала метрику r2. Конечно, это все отдельно для каждой прогнозирующей задачи.

Посмотрим, как признаки используются в моделях для каждого кластера. F-score на графиках ниже - частота появления признака в деревьях. Некоторые признаки почти не используются. Почистим от них. Threshold возьмем равным 0.01.

In [46]:
def cv_r2_score(n_splits, X_train, y_train, params={'objective': 'reg:squarederror', 'subsample': 1}):
    tscv = TimeSeriesSplit(n_splits=3)
    scores = []
    importances = []
    
    for tr_ix, cv_ix in tscv.split(X_train): 
        XGBR = myXGBR()
        XGBR.fit(X_train.iloc[tr_ix], y_train.iloc[tr_ix])
        y_pred = XGBR.predict(X_train.iloc[cv_ix])
        y_test = y_train.iloc[cv_ix].values
        scores.append(r2_score(y_test, y_pred))
        importances.append(XGBR.feature_importances_)
    score = np.mean(scores) # функция возвращает средне по cv значение r2
    feature_importances = np.array(importances).mean(axis=0) # и оценку важности признаков
    
    return score, feature_importances
In [47]:
%%time
features_dict = {}
dict_of_r2 = {}
for shift in range(1,7):
    data = shifted_data[shift]
    print('Модель T+{}; {}'.format(shift, datetime.now().strftime('%H:%M:%S')))
    features_dict[shift] = {}
    dict_of_r2[shift] = {}
    fig, ax = plt.subplots(1, 3, figsize=(15,5), dpi=500)
    for cluster_label, zone in enumerate(centroid_names):
        print('Кластер {}, зона {}; {}'.format(cluster_label, zone, datetime.now().strftime('%H:%M:%S')))
        # временной ряд региона
        timeseries = data.loc[:, [zone]]
        timeseries.columns = ['real_counts']
        timeseries.columns.name = zone
        # метаданные по региону
        #zone_data = marked_data[marked_data['PU_zone_id'] == zone].shift(-shift).dropna()
        
        # добавляем признаки
        timeseries_with_features = add_regressors(timeseries, K_s, K_h, K_d)
        timeseries_with_features.columns.name = zone
        
        # отделяем train
        timeseries_train = timeseries_with_features.loc[:split_date].copy()
        X_train = timeseries_train.drop(['real_counts'], axis=1)
        y_train = timeseries_train[['real_counts']]
        
        # оцениваем по кросс-валидации модель, обученную на всех признаках с параметрами по умолчанию
        r2_, feature_importances = cv_r2_score(4, X_train, y_train)
        dict_of_r2[shift][cluster_label] = [r2_]
        print(u'r2_score для всех признаков: {:.2f}'.format(r2_))
        
        # отбираем признаки из модели
        features_df = pd.DataFrame({'features': X_train.columns,
                                    'importances': feature_importances}).sort_values('importances',
                                                                                          ascending=False)

        selected_features = features_df[features_df['importances'] > 0.01]['features'].values
        features_dict[shift][cluster_label] = selected_features
        print('{} наиболее релевантных признаков'.format(len(selected_features)))

        # оцениваем по кросс-валидации модель, обученную на выбранных признаках с параметрами по умолчанию
        r2_, _ = cv_r2_score(4, X_train[selected_features], y_train)
        dict_of_r2[shift][cluster_label].append(r2_)
        print(u'r2_score для отобранных признаков: {:.2f}'.format(r2_))
        
        # визуализизуем релевантность отобранных признаков
        XGBR = myXGBR()
        XGBR.fit(X_train[selected_features], y_train)
        ax_ = ax[cluster_label]
        ax_ = plot_importance(XGBR, ax=ax_)
        ax_.set_title(u'Релевантность признаков для кластера {}'.format(cluster_label), size=12)
        ax_.set_ylabel('')    
        
        print('===')
        
    plt.subplots_adjust(wspace = .6, hspace = .3)
    plt.suptitle(u'Отбор признаков для модели T+{}'.format(shift))
    plt.show()
Модель T+1; 00:29:55
Кластер 0, зона 170; 00:29:55
r2_score для всех признаков: 0.95
10 наиболее релевантных признаков
r2_score для отобранных признаков: 0.94
===
Кластер 1, зона 79; 00:30:17
r2_score для всех признаков: 0.95
11 наиболее релевантных признаков
r2_score для отобранных признаков: 0.95
===
Кластер 2, зона 42; 00:30:41
r2_score для всех признаков: 0.74
11 наиболее релевантных признаков
r2_score для отобранных признаков: 0.68
===
Модель T+2; 00:31:07
Кластер 0, зона 170; 00:31:07
r2_score для всех признаков: 0.95
9 наиболее релевантных признаков
r2_score для отобранных признаков: 0.94
===
Кластер 1, зона 79; 00:31:32
r2_score для всех признаков: 0.95
13 наиболее релевантных признаков
r2_score для отобранных признаков: 0.95
===
Кластер 2, зона 42; 00:31:56
r2_score для всех признаков: 0.74
12 наиболее релевантных признаков
r2_score для отобранных признаков: 0.67
===
Модель T+3; 00:32:23
Кластер 0, зона 170; 00:32:23
r2_score для всех признаков: 0.95
7 наиболее релевантных признаков
r2_score для отобранных признаков: 0.94
===
Кластер 1, зона 79; 00:32:46
r2_score для всех признаков: 0.95
13 наиболее релевантных признаков
r2_score для отобранных признаков: 0.95
===
Кластер 2, зона 42; 00:33:11
r2_score для всех признаков: 0.74
10 наиболее релевантных признаков
r2_score для отобранных признаков: 0.69
===
Модель T+4; 00:33:37
Кластер 0, зона 170; 00:33:37
r2_score для всех признаков: 0.95
6 наиболее релевантных признаков
r2_score для отобранных признаков: 0.94
===
Кластер 1, зона 79; 00:33:57
r2_score для всех признаков: 0.95
11 наиболее релевантных признаков
r2_score для отобранных признаков: 0.95
===
Кластер 2, зона 42; 00:34:18
r2_score для всех признаков: 0.74
8 наиболее релевантных признаков
r2_score для отобранных признаков: 0.68
===
Модель T+5; 00:34:39
Кластер 0, зона 170; 00:34:39
r2_score для всех признаков: 0.95
8 наиболее релевантных признаков
r2_score для отобранных признаков: 0.94
===
Кластер 1, зона 79; 00:35:00
r2_score для всех признаков: 0.95
14 наиболее релевантных признаков
r2_score для отобранных признаков: 0.95
===
Кластер 2, зона 42; 00:35:23
r2_score для всех признаков: 0.74
11 наиболее релевантных признаков
r2_score для отобранных признаков: 0.69
===
Модель T+6; 00:35:46
Кластер 0, зона 170; 00:35:46
r2_score для всех признаков: 0.95
7 наиболее релевантных признаков
r2_score для отобранных признаков: 0.94
===
Кластер 1, зона 79; 00:36:07
r2_score для всех признаков: 0.95
12 наиболее релевантных признаков
r2_score для отобранных признаков: 0.95
===
Кластер 2, зона 42; 00:36:27
r2_score для всех признаков: 0.74
10 наиболее релевантных признаков
r2_score для отобранных признаков: 0.68
===
Wall time: 6min 53s

Даже если метрика не улучшилась, то уменьшение числа признаков ускоряет обучение.

In [48]:
with open('XGBoost/selected_features.pickle','wb') as file:
    pickle.dump(features_dict, file)
In [49]:
# посмотрим, как изменилась метрика r2 после отбора признаков
shift = 3
r2_df = pd.DataFrame(dict_of_r2[shift]).T
r2_df.columns=['all_features/default_params',
               'selected_features/default_params'
              ]
r2_df.index.name = 'cluster'

Например, так изменилась метрика R2 модели Т+3 при отборе признаков:

In [50]:
r2_df
# изменение метрики R2 модели Т+3 при отборе признаков
Out[50]:
all_features/default_params selected_features/default_params
cluster
0 0.945666 0.942118
1 0.949092 0.951052
2 0.735109 0.687903

6. Подбор гиперпараметров для XGBoost

Теперь надо подобрать гиперпараметры моделей каждого кластера. Будем оценивать модель с разными параметрами по метрике R2 на кросс-валидации через Bayesian Optimization Process (библиотека)

https://krasserm.github.io/2018/03/21/bayesian-optimization/ - разобраться при случае с этой статьей

In [51]:
from bayes_opt import BayesianOptimization
In [52]:
with open('XGBoost/selected_features.pickle', 'rb') as file:
    features_dict = pickle.load(file)
In [53]:
# функция, которую будем оптимизировать
def cv_score_opt(learning_rate, max_depth, gamma, n_estimators, reg_alpha, reg_lambda):
    xgb_params = {
        'objective': 'reg:squarederror',
        'subsample': 1,
        'learning_rate': learning_rate,
        'max_depth': int(max_depth),
        'gamma': gamma,
        'n_estimators': int(n_estimators),
        'reg_alpha': reg_alpha,
        'reg_lambda': reg_lambda
                 }
    # Cross-validation для временных рядов
    tscv = TimeSeriesSplit(n_splits=3)
    scores = []
    for tr_ix, cv_ix in tscv.split(X_train): 
        XGBR = myXGBR(**xgb_params)
        XGBR.fit(X_train.iloc[tr_ix], y_train.iloc[tr_ix])
        y_pred = XGBR.predict(X_train.iloc[cv_ix])
        y_test = y_train.iloc[cv_ix].values
        scores.append(r2_score(y_test, y_pred))
    score = np.mean(scores)
    
    return score
In [54]:
# сетка параметров
params_grid = {
    'learning_rate': (0,1),
    'max_depth': (1,100),
    'gamma': (0,10),
    'n_estimators': (10,300),
    'reg_alpha': (0,1),
    'reg_lambda': (0,1)
              }
In [55]:
%%time
best_hyperparams = {}
for shift in range(1,7):
    data = shifted_data[shift]
    print(f'Модель T+{shift}; {datetime.now().strftime("%H:%M:%S")}')
    best_hyperparams[shift] = {}
    for cluster_label, zone in enumerate(centroid_names):
        print(f'Кластер {cluster_label}, зона {zone}; {datetime.now().strftime("%H:%M:%S")}')
        
        # временной ряд региона
        timeseries = data.loc[:, [zone]]
        timeseries.columns = ['real_counts']
        timeseries.columns.name = zone
        
        # метаданные по региону - не используем
        #zone_data = marked_data[marked_data['PU_zone_id'] == zone].shift(-shift).dropna()
        
        # добавляем признаки
        timeseries_with_features = add_regressors(timeseries, K_s, K_h, K_d)
        timeseries_with_features.columns.name = zone
        
        # отделяем train и выбирыем нужные признаки
        selected_features = features_dict[shift][cluster_label]
        timeseries_train = timeseries_with_features.loc[:split_date].copy()
        X_train = timeseries_train.drop(['real_counts'], axis=1)[selected_features]
        y_train = timeseries_train[['real_counts']]

        # инициализируем BO оптимизатор
        optimizer = BayesianOptimization(cv_score_opt, params_grid, verbose=1)
        # максимизируем r2_score
        optimizer.maximize(n_iter=25, init_points=5, acq='ei')
        
        # подобранные лучшие параметры
        best_p = optimizer.max['params']
        best_p['max_depth'] = np.ceil(best_p['max_depth'])
        best_p['n_estimators'] = np.ceil(best_p['n_estimators'])
        
        # запишем в словарь
        xgb_params = {'objective': 'reg:squarederror',
                      'subsample': 1} # это константные параметры
        xgb_params.update(best_p)
        best_hyperparams[shift][cluster_label] = xgb_params
        print(xgb_params)

        # оцениваем по кросс-валидации модель, обученную на выбранных признаках с подобранными параметрами
        r2_, _ = cv_r2_score(4, X_train, y_train, params=xgb_params)
        print(u'r2_score для отобранных признаков c настроенными параметрами: {:.2f}'.format(r2_))
        dict_of_r2[shift][cluster_label].append(r2_)
        print('===')
Модель T+1; 00:36:49
Кластер 0, зона 170; 00:36:49
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  2        |  0.9379   |  9.252    |  0.129    |  33.65    |  99.22    |  0.1572   |  0.7365   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 9.252136390232458, 'learning_rate': 0.12900157849203775, 'max_depth': 34.0, 'n_estimators': 100.0, 'reg_alpha': 0.1571913311926394, 'reg_lambda': 0.7365036393427848}
r2_score для отобранных признаков c настроенными параметрами: 0.94
===
Кластер 1, зона 79; 00:44:54
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  4        |  0.9405   |  7.419    |  0.5385   |  82.37    |  263.3    |  0.2441   |  0.4059   |
|  5        |  0.9502   |  0.9248   |  0.2635   |  9.668    |  211.6    |  0.5265   |  0.2203   |
|  7        |  0.9502   |  2.749    |  0.3119   |  7.889    |  207.1    |  0.4663   |  0.7592   |
|  8        |  0.9548   |  7.455    |  0.07126  |  8.7      |  211.8    |  0.5103   |  0.0707   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 7.455408028683205, 'learning_rate': 0.07125831600885879, 'max_depth': 9.0, 'n_estimators': 212.0, 'reg_alpha': 0.5103226199012304, 'reg_lambda': 0.0707044949767981}
r2_score для отобранных признаков c настроенными параметрами: 0.95
===
Кластер 2, зона 42; 00:55:04
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  2        |  0.634    |  7.43     |  0.05845  |  60.43    |  222.7    |  0.03571  |  0.1525   |
|  3        |  0.6511   |  5.543    |  0.2472   |  47.84    |  238.2    |  0.7913   |  0.7433   |
|  15       |  0.6638   |  6.115    |  0.1585   |  47.15    |  238.6    |  0.6023   |  0.6636   |
|  16       |  0.668    |  6.575    |  0.1502   |  48.49    |  238.5    |  0.1391   |  0.6736   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 6.574660474124636, 'learning_rate': 0.15020221593020455, 'max_depth': 49.0, 'n_estimators': 239.0, 'reg_alpha': 0.1391496235163141, 'reg_lambda': 0.6736414973431143}
r2_score для отобранных признаков c настроенными параметрами: 0.68
===
Модель T+2; 01:08:34
Кластер 0, зона 170; 01:08:34
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  2        |  0.9294   |  1.498    |  0.3999   |  18.45    |  280.3    |  0.1477   |  0.1706   |
|  4        |  0.9323   |  0.6041   |  0.336    |  14.72    |  294.9    |  0.3351   |  0.6523   |
|  6        |  0.9434   |  9.516    |  0.2771   |  2.813    |  299.9    |  0.6404   |  0.4669   |
|  10       |  0.9436   |  7.481    |  0.09006  |  8.203    |  296.7    |  0.8091   |  0.3003   |
|  16       |  0.9437   |  2.624    |  0.2099   |  8.491    |  53.24    |  0.9592   |  0.5088   |
|  17       |  0.9472   |  5.378    |  0.1428   |  5.356    |  50.03    |  0.6965   |  0.7085   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 5.378336799220448, 'learning_rate': 0.14282064591859445, 'max_depth': 6.0, 'n_estimators': 51.0, 'reg_alpha': 0.6964803337538659, 'reg_lambda': 0.7085254848275877}
r2_score для отобранных признаков c настроенными параметрами: 0.94
===
Кластер 1, зона 79; 01:11:14
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 7.361199933903089, 'learning_rate': 0.11746905245053274, 'max_depth': 9.0, 'n_estimators': 148.0, 'reg_alpha': 0.19807496357042076, 'reg_lambda': 0.9153355059290995}
r2_score для отобранных признаков c настроенными параметрами: 0.95
===
Кластер 2, зона 42; 01:20:21
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  2        |  0.6343   |  3.416    |  0.0815   |  39.83    |  256.4    |  0.1385   |  0.2966   |
|  3        |  0.6556   |  1.344    |  0.1772   |  42.56    |  187.0    |  0.7207   |  0.8634   |
|  25       |  0.6591   |  1.504    |  0.04838  |  41.57    |  253.8    |  0.4401   |  0.7491   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 1.5037698052406812, 'learning_rate': 0.04838482927526366, 'max_depth': 42.0, 'n_estimators': 254.0, 'reg_alpha': 0.44012060906827033, 'reg_lambda': 0.7490924323557343}
r2_score для отобранных признаков c настроенными параметрами: 0.67
===
Модель T+3; 01:33:11
Кластер 0, зона 170; 01:33:11
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  2        |  0.9349   |  5.305    |  0.2378   |  13.36    |  229.7    |  0.8655   |  0.709    |
|  28       |  0.9446   |  3.288    |  0.2349   |  3.481    |  230.3    |  0.4623   |  0.6584   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 3.2880318519599125, 'learning_rate': 0.23489195953757924, 'max_depth': 4.0, 'n_estimators': 231.0, 'reg_alpha': 0.4622855004640861, 'reg_lambda': 0.6584319303076308}
r2_score для отобранных признаков c настроенными параметрами: 0.94
===
Кластер 1, зона 79; 01:37:28
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  4        |  0.9463   |  7.638    |  0.3944   |  85.32    |  124.8    |  0.8221   |  0.5073   |
|  8        |  0.9492   |  9.826    |  0.13     |  84.61    |  122.7    |  0.4777   |  0.9149   |
|  16       |  0.9498   |  6.995    |  0.05422  |  86.1     |  123.4    |  0.3425   |  0.6736   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 6.995497175896174, 'learning_rate': 0.05421684938179916, 'max_depth': 87.0, 'n_estimators': 124.0, 'reg_alpha': 0.34247173735774705, 'reg_lambda': 0.6735705522889324}
r2_score для отобранных признаков c настроенными параметрами: 0.95
===
Кластер 2, зона 42; 01:48:01
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  3        |  0.6561   |  6.228    |  0.09421  |  75.77    |  132.1    |  0.8198   |  0.2114   |
|  9        |  0.6569   |  6.12     |  0.1793   |  71.32    |  133.7    |  0.4287   |  0.2486   |
|  13       |  0.7159   |  7.764    |  0.0702   |  5.563    |  133.6    |  0.8498   |  0.8426   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 7.763976196760593, 'learning_rate': 0.0701982562695197, 'max_depth': 6.0, 'n_estimators': 134.0, 'reg_alpha': 0.8498392865069826, 'reg_lambda': 0.8425766259155448}
r2_score для отобранных признаков c настроенными параметрами: 0.69
===
Модель T+4; 01:55:18
Кластер 0, зона 170; 01:55:18
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  2        |  0.9326   |  8.735    |  0.2695   |  46.67    |  262.9    |  0.3088   |  0.8423   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 8.734929616344132, 'learning_rate': 0.26949403169819885, 'max_depth': 47.0, 'n_estimators': 263.0, 'reg_alpha': 0.3088387318987521, 'reg_lambda': 0.8423476751668452}
r2_score для отобранных признаков c настроенными параметрами: 0.94
===
Кластер 1, зона 79; 02:01:58
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  2        |  0.9482   |  2.21     |  0.04982  |  54.44    |  91.29    |  0.9866   |  0.3164   |
|  20       |  0.9486   |  0.397    |  0.04099  |  59.33    |  165.6    |  0.8392   |  0.6853   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 0.3970175241127938, 'learning_rate': 0.040990470606331986, 'max_depth': 60.0, 'n_estimators': 166.0, 'reg_alpha': 0.8392366188123002, 'reg_lambda': 0.6853006414982076}
r2_score для отобранных признаков c настроенными параметрами: 0.95
===
Кластер 2, зона 42; 02:10:40
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  2        |  0.6425   |  8.857    |  0.09295  |  88.47    |  104.0    |  0.3603   |  0.4696   |
|  16       |  0.6943   |  5.181    |  0.6063   |  2.308    |  208.9    |  0.9224   |  0.4315   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 5.180562947819044, 'learning_rate': 0.6062669355057018, 'max_depth': 3.0, 'n_estimators': 209.0, 'reg_alpha': 0.9223508573356712, 'reg_lambda': 0.4315031483789883}
r2_score для отобранных признаков c настроенными параметрами: 0.68
===
Модель T+5; 02:18:39
Кластер 0, зона 170; 02:18:39
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  2        |  0.93     |  2.783    |  0.3046   |  37.76    |  125.6    |  0.8763   |  0.112    |
|  5        |  0.9359   |  9.256    |  0.2071   |  89.59    |  61.93    |  0.9628   |  0.9387   |
|  12       |  0.9374   |  6.029    |  0.1503   |  38.09    |  229.8    |  0.2033   |  0.6284   |
|  29       |  0.9387   |  5.776    |  0.06424  |  39.33    |  232.3    |  0.005905 |  0.6313   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 5.7763573297790325, 'learning_rate': 0.06423983067546557, 'max_depth': 40.0, 'n_estimators': 233.0, 'reg_alpha': 0.005905207386423528, 'reg_lambda': 0.6312731526973356}
r2_score для отобранных признаков c настроенными параметрами: 0.94
===
Кластер 1, зона 79; 02:25:48
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  2        |  0.95     |  8.095    |  0.1095   |  31.6     |  175.1    |  0.9436   |  0.3861   |
|  3        |  0.9505   |  1.636    |  0.116    |  46.2     |  165.7    |  0.0214   |  0.3705   |
|  13       |  0.9518   |  7.434    |  0.09353  |  18.98    |  65.06    |  0.2015   |  0.9674   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 7.433558668833672, 'learning_rate': 0.0935278117297772, 'max_depth': 19.0, 'n_estimators': 66.0, 'reg_alpha': 0.20151513506205532, 'reg_lambda': 0.967390572907068}
r2_score для отобранных признаков c настроенными параметрами: 0.95
===
Кластер 2, зона 42; 02:31:24
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  2        |  0.5008   |  2.288    |  0.951    |  43.14    |  79.97    |  0.6483   |  0.5577   |
|  3        |  0.5229   |  7.968    |  0.8944   |  8.283    |  224.7    |  0.2798   |  0.3313   |
|  4        |  0.6506   |  7.465    |  0.2954   |  67.58    |  143.2    |  0.1621   |  0.5969   |
|  8        |  0.6787   |  8.309    |  0.05787  |  36.58    |  264.8    |  0.3808   |  0.9277   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 8.308948394038278, 'learning_rate': 0.05787253644236612, 'max_depth': 37.0, 'n_estimators': 265.0, 'reg_alpha': 0.38079576837869, 'reg_lambda': 0.927668941884226}
r2_score для отобранных признаков c настроенными параметрами: 0.69
===
Модель T+6; 02:40:00
Кластер 0, зона 170; 02:40:00
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  4        |  0.9345   |  6.3      |  0.2356   |  39.07    |  291.8    |  0.0612   |  0.6636   |
|  11       |  0.9368   |  6.76     |  0.1063   |  67.79    |  147.6    |  0.924    |  0.4953   |
|  19       |  0.9383   |  0.2761   |  0.06477  |  13.35    |  48.06    |  0.2903   |  0.8618   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 0.276124329784625, 'learning_rate': 0.06476970272843974, 'max_depth': 14.0, 'n_estimators': 49.0, 'reg_alpha': 0.2903306760606683, 'reg_lambda': 0.861822018476941}
r2_score для отобранных признаков c настроенными параметрами: 0.94
===
Кластер 1, зона 79; 02:45:04
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  3        |  0.9523   |  6.558    |  0.6051   |  2.672    |  163.6    |  0.6803   |  0.4717   |
|  16       |  0.9537   |  3.868    |  0.5649   |  2.449    |  272.1    |  0.9917   |  0.6128   |
|  21       |  0.9545   |  2.639    |  0.1688   |  3.081    |  271.9    |  0.6848   |  0.8067   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 2.6393937547883173, 'learning_rate': 0.1688479055842439, 'max_depth': 4.0, 'n_estimators': 272.0, 'reg_alpha': 0.6847711549656149, 'reg_lambda': 0.8066752873135591}
r2_score для отобранных признаков c настроенными параметрами: 0.95
===
Кластер 2, зона 42; 02:49:37
|   iter    |  target   |   gamma   | learni... | max_depth | n_esti... | reg_alpha | reg_la... |
-------------------------------------------------------------------------------------------------
|  2        |  0.6419   |  1.888    |  0.3346   |  30.52    |  145.6    |  0.4043   |  0.3856   |
|  8        |  0.6989   |  2.092    |  0.05631  |  8.159    |  198.9    |  0.003237 |  0.5757   |
|  26       |  0.7035   |  6.16     |  0.2028   |  4.414    |  251.7    |  0.5881   |  0.9466   |
=================================================================================================
{'objective': 'reg:squarederror', 'subsample': 1, 'gamma': 6.160471735286376, 'learning_rate': 0.20282845500869717, 'max_depth': 5.0, 'n_estimators': 252.0, 'reg_alpha': 0.5881201364183917, 'reg_lambda': 0.9466382371793499}
r2_score для отобранных признаков c настроенными параметрами: 0.68
===
Wall time: 2h 19min 6s
In [56]:
with open('XGBoost/best_hyperparams.pickle','wb') as file:
    pickle.dump(best_hyperparams, file)
In [57]:
# посмотрим, как изменилась метрика r2 после настройки гиперпараметров модели
shift = 3
r2_df = pd.DataFrame(dict_of_r2[shift]).T
r2_df.columns=['all_features/default_params',
               'selected_features/default_params',
               'selected_features/optimized_params'
              ]
r2_df.index.name = 'cluster'
In [58]:
r2_df
Out[58]:
all_features/default_params selected_features/default_params selected_features/optimized_params
cluster
0 0.945666 0.942118 0.942118
1 0.949092 0.951052 0.951052
2 0.735109 0.687903 0.687903

7. Прогнозирование на июль 2018 года

Выбранными моделями построим для каждой географической зоны и каждого конца истории от 2018.06.30 23:00 до 2018.07.31 17:00 прогнозы на 6 часов вперёд; посчитаем ноутбуке ошибку прогноза по следующему функционалу:

$$Q_{july2018}=\frac{1}{R∗739∗6}\sum_{r=1}^{R}\sum_{T=2018.06.30 23:00}^{2018.07.31 17:00}\sum_{i=1}^{6}∣\widehat{y}^{r}_{T|T+i}−{y}^{r}_{T+i}| $$

, где R - это число рядов (83 у нас), 739 - количество часов в прогнозируемом периоде.

In [59]:
# прогнозируем на июнь
july = pd.date_range(datetime(2018,6,30,23), datetime(2018,7,31,23), freq='H')

hour = pd.DateOffset(hours=1)
In [60]:
len(july)-6
Out[60]:
739
In [61]:
with open('XGBoost/selected_features.pickle', 'rb') as file:
    features_dict = pickle.load(file)
with open('XGBoost/best_hyperparams.pickle', 'rb') as file:
    hyperparams_dict = pickle.load(file)
In [62]:
%%time
july_forecasts = []
for shift in range(1,7):
    data = shifted_data[shift]
    print(f'Модель T+{shift}; {datetime.now().strftime("%H:%M:%S")}')
    best_hyperparams[shift] = {}
    for zone in zone_cluster_labels.keys():
        cluster_label = zone_cluster_labels[zone]
        #print(f'Зона {zone}; {datetime.now().strftime("%H:%M:%S")}')
        
        # временной ряд региона
        timeseries = data.loc[:, [zone]]
        timeseries.columns = ['real_counts']
        timeseries.columns.name = zone
        
        # метаданные по региону
        #zone_data = marked_data[marked_data['PU_zone_id'] == zone].shift(-shift).dropna()
        
        # добавляем признаки
        timeseries_with_features = add_regressors(timeseries, K_s, K_h, K_d)
        timeseries_with_features.columns.name = zone
        
        # нужные признаки
        selected_features = features_dict[shift][cluster_label]
        
        ## подобранные гиперпараметры модели
        xgb_params = hyperparams_dict[shift][cluster_label]
        xgb_params.update({'max_depth': int(xgb_params['max_depth']),
                           'n_estimators': int(xgb_params['n_estimators'])})

        #### ПРОГНОЗИРУЕМ НА ИЮЛЬ
        ## отделяем train и test
        timeseries_train = timeseries_with_features.loc[:split_date].copy()
        X_train = timeseries_train.drop(['real_counts'], axis=1)[selected_features]
        y_train = timeseries_train[['real_counts']]
        
        timeseries_test = timeseries_with_features.loc[july[0]:july[-7]].copy()
        X_test = timeseries_test.drop(['real_counts'], axis=1)[selected_features]
        y_test = timeseries_test[['real_counts']]  
        
        # модель
        XGBR = myXGBR(**xgb_params)
        XGBR.fit(X_train, y_train)
        ## прогноз
        july_forecast = pd.DataFrame(XGBR.predict(X_test).ravel(),
                                     index = X_test.index,
                                     columns=['prediction'])
        july_forecast = pd.concat([y_test, july_forecast], axis=1)
        july_forecast['zone'] = zone
        july_forecast['shift'] = shift
        july_forecasts.append(july_forecast)
Модель T+1; 02:55:56
Модель T+2; 03:08:00
Модель T+3; 03:15:45
Модель T+4; 03:20:55
Модель T+5; 03:35:00
Модель T+6; 03:53:22
Wall time: 59min 39s
In [63]:
july_forecasts = pd.concat(july_forecasts, axis = 0)
In [64]:
july_forecasts.to_csv('XGBoost/july_forecasts.csv', index=True)
In [65]:
# посчитаем Q_july
Q_july = np.sum(np.abs(july_forecasts['real_counts'] - july_forecasts['prediction']))/ (R*739*6)
print('Q_july= {:.2f}'.format(Q_july)) 
Q_july= 15.43

Кажется, это не очень большая средняя ошибка прогноза в масштабах спроса на такси в Нью-Йорке.